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Abstract 

In this work, numerical solutions of the two dimensional time dependent incompress- 
ible flow, in a driven cavity at high Reynolds number Re, are presented. At high 
Re, there is a controversy. Some studies predicted that the flow is steady, others 
found time dependent non-steady flow, either periodic or aperiodic. In this study, 
the driven lid cavity is successfully solved using a very fine grid mesh, for Re up to 
30 000. We discretize the Vorticity-Stream formulation of the Navier-Stokes equa- 
tion with the SSPRK(5,4) scheme in a 1024 x 1024 grid. Using this very fine grid, 
the results obtained converge to a stationary solution. Detailed results for Re be- 
tween 5 000 and 30 000 are presented. The driven lid cavity problem is solved with 
a NVIDIA GPU using the CUDA programming environment with double precision. 

Key words: Driven Lid Cavity Flow, 2-D Time Dependent Incompressible N-S 
Equation, Fine Grid, Reynolds Number, CUDA 



1 Introduction 

The driven lid cavity flow is one of the most studied problems of the physics of 
fluids. The simple geometry makes the problem easy to solve numerically and 
apply boundary conditions. Despite being a problem rather studied, there are 
still some questions and controversy about what happens at high Reynolds 
numbers. 

In the literature it is possible to flnd numerous studies about the driven cavity 
flow, however the nature of the flow at high Reynolds number is still not agreed 
upon. Driven cavity flow serve as a benchmark problem for numerical meth- 
ods in terms of numerical efficiency and accuracy. Erturk [1] grouped these 
numerical studies into three categories. In the first category are the studies 
with numerical solutions of 2-D steady incompressible fiow at high Reynolds 
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numbers. Some of these studies are Erturk et. al. [2], Erturk and Gokcol [3], 
Barragy and Carey [4], Schreiber and Keller [5], Benjamin and Denny [6], Liao 
and Zhu [7], Ghia et. al. [8] for Re = 10 000. Barragy and Carey [4] have also 
presented solutions for Re = 12 500. Erturk et. al. [2] and also Erturk and 
Gokcol [3] have presented steady solutions up to Re = 30 000. In the second 
category we have the studies about hydrodynamic stability analysis, Fortin 
et. al. [9], Gervais et. al. [10], Sahin and Owens [11] and Abouhamza and 
Pierre [12]. And the third category includes the studies about the steady to 
unsteady transition flow through a Direct Numerical Simulation (DNS) and 
the transition Reynolds number. Some of the studies found in the literature are 
Auteri, Parolini and Quartepelle [13], Peng, Shiau and Hwang [14], Tiesinga, 
Wubs and Veldman [15], Poliashenko and Aidun [16], Cazemier, Verstappen 
and Veldman [17], Goyon [18], Wan, Zhou and Wei [19] and LiflFman [20]. 

Erturk [21] studied the numerical solutions of 2-D steady flow, with time 
independent equations, for Re < 20 000 and found that at high Reynolds 
numbers the solutions obtained depend of the mesh spacing. For coarse meshes, 
at high Reynolds numbers, he found no convergence. But for a flne grid he 
flnds a steady solution. Thus he claimed that the true solution is steady and 
to be found numerically a very flne grid is needed. 

Zhen-Hua et al. [22] with multi-relaxation-time (MRT model) lattice Boltz- 
mann method obtained solutions for up to Re = 1 000 000, but the grid used 
(256 X 256) in the MRT model are coarse compared to the article of Erturk, 
and the solutions obtained by them are not steady. 

Notice that the existence of solutions to the stationary flow equation, does not 
necessarily imply that the natural solution is stationary. Let us illustrate this 
claim with the well known case of metastable water. At the normal pressure 
and at temperatures below freezing, two phases of water, i.e., two solutions of 
the non-linear equations for the water, exist. The stable solution is composed 
of solid water (ice) and the metastable is composed of supercooled water. In 
the same way the existence of a stationary solution to the flow equations, 
does not necessarily exclude the existence of a second, unsteady, solution. The 
solution of the time-dependent equation is necessary to discriminate which is 
the true solution, occurring naturally in the lab. 

Since diflFerent authors flnd diflFerent solutions, in this paper we examine if for 
high Reynolds numbers as high as Re = 30 000 the numerical solutions of 2-D 
time dependent incompressible flow in a driven cavity become steady after a 
given time. Importantly, our equations are time-dependent, and our grid is 
very flne, thus we expect to obtain the true solution. We use the vorticity- 
stream function formulation to the Navier-Stokes equation. 

This study is divided in flve sections. The flrst section, the objective is pre- 
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sented and the existing works on the subject. In the second section, we present 
the theory and the fundamental equations used in the numeric calculations. 
In the third section, we describe the driven lid cavity problem, numeric equa- 
tions, boundary conditions and the algorithm used. In the section four, we 
present the results and a comparison with existing studies about the cavity 
flow. Finally, in the section flve, we present the conclusions about this study. 



2 Vorticity-Stream Function Formulation 



The fundamental equations of fluid dynamics are based on three universal 
laws of conservation: mass, moment and energy conservation. In general, there 
are two ways to solve numerically the Navier-Stokes equations. The flrst is 
proceeding with primitive variables formulation, v e p. The second is using 
the vorticity-stream function approach. 

Using the stream function, i/j and the vorticity, uu^ in place of the primitive 
variables v e where these quantities are deflned by: 

^ ^ dv du 

and 

I (2) 

Combining the equations (1) and (2), 

A^ = V^=— + — = -c. (3) 



The main reason for enter the stream function is that, for runoflFs in which p 
and V are constants, the equation of continuity is satisfled. 



The vorticity-stream function is given by 

duo duo dip duo 

dt dy dx dx dy 
d'^ip d'^ip 
dx'^ dy'^ 

where ~u = (ix, v) is the fluid velocity, p the pressure, p the fluid density and 



(4) 
(5) 
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V kinematic viscosity. The kinematic viscosity is given by: 

r] 



P 



(6) 



where t] is the viscosity. 



3 Numerical Method 



In this section we present the driven hd cavity problem, the numerical equa- 
tions, the boundary conditions and the algorithm. In the figure 1 is repre- 
sented the outline problem to study. The top of the wall cavity moves with 
speed u = U. 




Figure 1. Driven Lid Cavity. 



Numeric vorticity equation is given by: 



with 



(7) 



L K) = 



CO, 



-u. 



/ 1 » t' 

2Ax 



2Ay 



{Axf 



(Jj, 



(8) 



The five stage fourth order SSPRK developed by Ruuth and Spiteri, [23] , and 
guaranteed optimal [24], is given by: 
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uj(^) = cu" + 0.39175222657 At L (cu") (9) 
J2) ^ 0.444370493651235 cu" + 0.555626506348765 u^^^ (10) 
+ 0.368410593050371 At L (u;^^)) 

= 0.620101851488403 a;" + 0.379898148511597 u;^^) (11) 
+ 0.251891774271694 At L (u^^^) 

= 0.178079954393132 o;'^ + 0.821920045606868 a;^^) (12) 
+ 0.544974750228521 At L (u^^^) 
a;"+i = 0.517231671970585 (13) 
+ 0.096059710526147 + 0.063692468666290 AtL (cj^^^) 
+ 0.386708617503269 + 0.226007483236906 AtL (ic;^^)) 



To get the stream function, we apply the fourth order method, [3], to solve 
the Poisson equation (3), 



To solve the fourth order Stream equation, we use the Successive Over Re- 
laxation method (SOR), since this method converge much faster than the 
traditional methods (methods of Jacobi and Gauss-Seidel). The idea is to use 
the available current estimates from other locations (V'l+ij) when they be- 
come available and use the estimates of current positions (i^ij) to improve the 
estimative of ipfj^'- 



^n+i = ^(l^A + B + C) + il-p) ijl^ (15) 

O Qj 



where (3 is the relaxation parameter, and 
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^ = b (C+i, + C-^li) + c {^P:J+, + ^P:Jt\) (16) 

B = ^ + u^^^, + Su^ + + u:l_,) (17) 

C = (^IVi.+i - 2^+1 + ^ri+i - 2^IVi.- (18) 

a = 6 + c (19) 



To obtain the velocity {u e v) it is only need to numerically solve the equation 
(2), 



u. 



2Ay 
2Ax 



(22) 



3.1 Boundary Conditions 



The boundary conditions in the four walls are given by: 
• left and right walls: 

^ 'q^ Uyjall 



s = o 



UJ — 



• top and bottom walls: 



^ ~dy ^wall 

v = 

- = -0 



(23) 



(24) 



where U^aii is the uniform velocity for the translation of the top wall and zero 
for the other three walls. 
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In order to solve equation (13), Thorn's method is used for calculating vorticity 
on the boundaries, therefore the boundary conditions are given by 



< 

3.2 Algorithm 

An algorithm for solving the problem of the driven lid cavity using the vorticity- 
stream function approach is given by the scheme in figure 2. 
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Ax2 



A2/2 



^Ax 



(25) 



I nitial 
Conditions 



■ If not converge 



Solve Poisson 
Equation 



Ol3tain Velocity 
u and V 



Solve Vorticity 
Equation 



If converge 



See Results 



Figure 2. Scheme for solving the problem of the driven lid cavity using the vortici- 
ty-stream function approach. 



To ensure the stability and convergence of the algorithm. At should be small 
enough to a given viscosity, z/, and resolution of the grid. Ax Ay. The Reynolds 
number. Re, can be calculated using the kinematic viscosity, z/, and the con- 
ditions of the cavity: 

UwallLy pU^aiiLy 

L = LxLy 
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in the following calculations it was considered: 



P 
L 




1 



'y 



Using this values, it's the same using dimensionless variables. 



In this study, we defined RES by the following equations 



RESn 



RES, 



- y k'^^ 

N ^ I ^''.J 





(26) 



(27) 



for stream function and vorticity respectively, as convergence criteria to the 
steady state. 



4 Results 

In this work, we present results of the cavity fiow from Re = 5 000 up to 
Re = 30 000. 

The numerical code was written for use in CPU's (Intel(R) Core(TM)2 Quad 
CPU Q9450 @ 2.66GHz) and in CPU's (NVIDIA CEFORCE 280 CTX CPU 
with double precision capabilities). For the CPU, we use OPENMP with C 
language and for CPU we use CUDA language. Most of the CPU's only sup- 
port single precision but the recent CPU's have included support for double 
precision, although using double precision is almost eight times slower than 
single precision. Nevertheless, the CPU used is faster than the most recent 
Intel Quad core. Figure 3 summarizes our CPU code performance relative to 
the serial CPU version of our code and CPU parallel code performance relative 
to the serial CPU version of our code using double precision in all cases. In the 
parallel CPU version code we tested it in a Quad Core CPU with OPENMP. 

In order to obtain a good convergence of the method we need to choose an 
appropriate value for the (5 parameter, the relaxation parameter, in equation 
(15). We make several tests in order to see which value is the best and this 
parameter varies from 0.8 to 0.1 as the Reynolds number increase. 

Using the algorithm described above, we started with a grid of 220 x 220 points 
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Figure 3. GPU code speedup relative to the serial and parallel CPU code for a grid 
of 1024x1024. 



Re 


Grid 


At (s) 


Figures 


5 000 


1024 X 1024 


0.001 


5 and 11a 


10 000 


1024 X 1024 


0.001 


6 and lib 


15 000 


1024 X 1024 


0.001 


7 and 11c 


20 000 


1024 X 1024 


0.001 


8 and lid 


25 000 


1024 X 1024 


0.001 


9 and lie 


30 000 


1024 X 1024 


0.001 


10 and llf 



Table 1 
Grid mesh. 

for Re = 5 000 and Re = 10 000 and the solution converges to a steady state 
in time. 

For Re = 15 000 we obtain a solution that converges to a steady state with a 
grid of 300 x 300. With this small lattice we cannot obtain a solution for the 
Reynolds number above 15 000. Erturk et. al. [2] said that they have to use a 
grid of 1025 x 1025 for Re > 15 000. So, when we use a grid of 1024 x 1024 
for Re > 15 000, we obtain a steady solution up to Re = 30 000. The choice 
of this grid is due to the GPU architecture. Table 1 shows the grid mesh used 
at various Reynolds numbers for the final results. For Re up to 20 000, the 
algorithm ended when the value obtained by RES^ and RES^ were less than 
10~^^ and 10~^^ for stream function and vorticity respectively. For Re = 25 000 
and Re = 30 000, the values considered for RES^ and RES^ were 10~^° and 
10~^. Such values are more than satisfactory, demonstrating that the solution 
converges to a steady state. 

The figure 4 shows the boundary conditions and a schematic of the vortices 
generated in the driven cavity fiow. In this figure, the abbreviations TL, BL 
and BR refer to top left, bottom left and bottom right corners of the cavity, 
respectively, and the number following these abbreviations to the vortices that 
appear in the fiow, numbered according to size, Erturk [21]. 
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BL2\ ^ 




^ BR2 








/ / 



u=0 v=0 



Figure 4. Schematic view of the driven hd cavity flow. Source: Erturk [21]. 



Figures 5 to 9a show the stream function contours and Figure 11 the vorticity 
contours of the cavity flow up to Re = 30 000. In the flgure 5 are represented 
the stream function contours for Re = 5 000, according to flgure 4, all the 
vortices appear except the BL3 and TL2. This result is in agreement with the 
result obtained by Erturk [1]. For Re = 10 000, flgure 6, all the vortices appear 
except the TL2. For Re > 15 000, all the vortices can be seen. These contour 
flgures show that, the flne grid mesh provides very smooth solutions at high 
Reynolds numbers. 
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0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 



0.955 0.96 0.965 0.97 0.975 0.98 0.985 0.99 0.995 1 



(b) (c) 

Figure 5. (a) Stream function contours for Re = 5 000 with 1024 x 1024 points, (b) 
and (c) correspond to the Stream function in the bottom left and right corners of 
the cavity respectively. The stream function is expressed in m?s~^ and the sides of 
the cavity in m. 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

(a) 




0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 0.95 0.955 0.96 0.965 0.97 0.975 0.98 0.985 0.99 0.995 1 



(b) (c) 

Figure 6. (a) Stream function contours for Re = 10 000 with 1024 x 1024 points, (b) 
and (c) correspond to the Stream function in the bottom left and right corners of 
the cavity respectively. The stream function is expressed in m'^s~^ and the sides of 
the cavity in m. 
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Figure 7. (a) Stream function contours for Re = 15 000 with 1024 x 1024 points, (b) 
and (c) correspond to the Stream function in the bottom left and right corners of 
the cavity respectively. The stream function is expressed in im?s~^ and the sides of 
the cavity in m. 
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-5.5E-005 
-6E-005 
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-7E-005 
-7.5E-005 



0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 



(b) 




0.95 0.955 0.96 0.965 0.97 0.975 0.98 0.985 0.99 0.995 1 



(c) 



Figure 8. (a) Stream function contours for Re = 20 000 with 1024 x 1024 points, (b) 
and (c) correspond to the Stream function in the bottom left and right corners of 
the cavity respectively. The stream function is expressed in im?s~^ and the sides of 
the cavity in m. 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

(a) 




0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 0.95 0.955 0.96 0.965 0.97 0.975 0.98 0.985 0.99 0.995 1 



(b) (c) 

Figure 9. (a) Stream function contours for Re = 25 000 with 1024 x 1024 points, (b) 
and (c) correspond to the Stream function in the bottom left and right corners of 
the cavity respectively. The stream function is expressed in m?s~^ and the sides of 
the cavity in m. 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

(a) 




0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05 0.95 0.955 0.96 0.965 0.97 0.975 0.98 0.985 0.99 0.995 1 



(b) (c) 

Figure 10. (a) Stream function contours for Re = 30000 with 1024 x 1024 points, 
(b) and (c) correspond to the Stream function in the bottom left and right corners 
of the cavity respectively. The stream function is expressed in m?s~^ and the sides 
of the cavity in m. 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



(a) Re = 5000 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

(b) Re = 10 000 




5 Conclusion 



Numerical solutions of 2-D time dependent incompressible flow at high Reynolds 
numbers in a driven cavity are presented. The numerical equations are solved 
computationally using the numerical method above. For the Reynolds num- 
bers studied, up to 30 000, the solution converges to a stationary state when 
using a very flne grid. 

Based on the papers of Erturk et. al. [2], Erturk [1] and in this study, we 
conclude that in order to obtain a steady solution for the driven cavity flow, a 
very flne grid mesh is necessary when high Reynolds numbers are considered 
and also at high Reynolds numbers when a coarse grid mesh is used then the 
solution oscillates. This happens because a coarse mesh is not able to include 
the very small vortices at the corners. 

According to Erturk [1], the studies that presented unsteady solutions of driven 
cavity flow using Direct Numerical Simulations (DNS) ([13,17,18,20,14,16,15,19]) 
have experienced the same type of numerical oscillations because they have 
used a coarse mesh. With this study we agreed with Erturk [1]. In all of 
the Direct Numerical Simulation studies on the driven cavity flow found 
in the literature ([13,17,18,20,14,16,15,19]), the maximum number of grid 
points used is less than 300 x 300. Therefore, the periodic solutions found 
in [13,17,18,20,14,16,15,19] are similar to the false periodic solutions observed 
by Erturk et. al. [2] and Erturk [1] when a coarse grid mesh is used. 

In short, if a suflBciently flne grid mesh is used, a Direct Numerical Simulation 
algorithm or other algorithm would give the same steady results obtained by 
Erturk et. al. [2], Erturk [1] and in this study. 
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